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Abstract 

In this paper, we propose a method for the eonstruetion of loeally eonservative flux fields through 
a variation of the Generalized Multiseale Finite Element Method (GMsFEM). The flux values are 
obtained through the use of a Ritz formulation in whieh we augment the resulting linear system of 
the eontinuous Galerkin (CG) formulation in the higher-order GMsFEM approximation spaee. In 
partieular, we impose the finite volume-based restrietions through ineorporating a sealar Eagrange 
multiplier for eaeh mass eonservation eonstraint. This approaeh ean be equivalently viewed as 
a eonstraint minimization problem where we minimize the energy funetional of the equation re- 
strieted to the subspaee of funetions that satisfy the desired eonservation properties. To test the 
performanee of the method we eonsider equations with heterogeneous permeability eoeffieients 
that have high-variation and diseontinuities, and eouple the resulting fluxes to a two-phase flow 
model. The inerease in aeeuraey assoeiated with the eomputation of the GMsFEM pressure solu¬ 
tions is inherited by the flux fields and saturation solutions, and is elosely eorrelated to the size of 
the redueed-order systems. In partieular, the addition of more basis funetions to the enriehed mul¬ 
tiseale spaee produees solutions that more aeeurately eapture the behavior of the fine seale model. 
A variety of numerieal examples are offered to validate the performanee of the method. 
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1. Introduction and problem statement 

In this paper, we primarily eonsider the equation given by 


— div{Ak{x)'Vp) = q in 

P = Pd on r^, (1) 

—AkVp-n = gN on F^r 


where k{x) is a heterogeneous field with high contrast. In particular, we assume that there is 
a positive constant fcmin such that k{x) > fcmin > 0, while k{x) can have very large values (i.e., 
fcmax /krain IS large). Additionally, A is a known mobility coefficient, q denotes any external forcing, 
and p is an unknown pressure field satisfying Dirichlet and Neumann boundary conditions given by 
Pd and gj^, respectively. Here kl a convex polygonal and two dimensional domain with boundary 
dkl = F £) U F jv. 

Let us consider a function in H^{D) whose trace on F^) coincides with the given value we 
denote this function also by p^. The variational formulation of problem O is to find p G 
with (p — pd) G = {tu G = 0} and such that 

a(p, v) = F{v) - {gN, 'u)r^ for all v G H^, (2) 


where, for p, u G the bilinear form a is defined by 


a{p,v) = / Ak{x)Vp{x)Vv{x)dx, 


the functional F is defined by 

F{v) = / q{x)v{x)dx 

Jn 

and the linear functional related to the boundary condition is given by 


( 3 ) 

( 4 ) 


{gN,v)rN 


gj^{x)v{x)dl. 


( 5 ) 


A main goal of our work is to obtain conservative discretizations of the equations above. More 
specifically, construction of approximation that satisfy some given conservation of mass restriction 
on subdomains of interest. We note that a popular conservative discretization is the Finite Volume 
(FV) method. The classical FV discretization provides and approximation of the solution in the 
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space of piecewise linear functions with respect to a triangulation while satisfying conservation of 
mass on elements of a dual triangualtion. When the approximation of the piecewise linear space 
is not enough for the problem at hand, advance approximation spaces need to be used. However, 
in some cases this requires a sacrifice of the conservation properties of the FV method. In this 
work we present an extension of the FV method for general approximation spaces that enrich 
classical approximation spaces (such as the space of piecewise linear functions). In particular, this 
conservative discretization can be used in conjunction with recently introduced GMsFEM spaces. 

We note that FV methods that use higher degree piecewise polynomials have been introduced 
in the literature. The fact that the dimension of the approximation spaces is larger than the number 
of restrictions led the researchers of [l], Q] to introduce additional control volumes to match the 
number of restrictions to the number of unknowns. An alternative approach is to consider a Petrov- 
Galerkin formulation with additional test functions rather that only piecewise constant functions 
on the dual grid. They where able to obtain stability of the method as well as error estimates. It 
is important to observe that the additional control volumes require additional computational work 
to be constructed and in some cases are not easy to construct (see also [3, ^). Additionally, it 
is well know that piecewise smooth approximations spaces do not perform well for multiscale 
high-contrast problems ll5Ul2l]. Another technique that one can use in order obtain a conservative 
method with richer approximations spaces is the following (see for instance OOll where the authors 
use a similar approach). In the discrete linear system obtained by a finite element discretization, 
it is possible to substitute appropriate number of equations by finite volume equations involving 
only the standard dual grid. This approach has the advantage that no additional control volume 
needs to be constructed. It may have the flexibility of both FV and FE procedures given as mass 
conservative fluxes and residual minimization properties. Some previous numerical experiments 
suggest a drawback of this approach - the resulting discrete problems may be ill-conditioned for 
large dimension coarse spaces, specially for higher order approximation spaces and multiscale 
problems. 

In this paper we propose the alternative of using a Ritz formulation and construct a solution 
procedure that combines a continuous Galerkin-type formulation that concurrently satisfies mass 
conservation restrictions. To this end, we augment the resulting linear system of the Galerkin 
formulation in the higher order approximation space to impose the finite volume restrictions. In 
particular, we do that by using a scalar Lagrange multiplier for each restriction. This approach 
can be equivalently viewed as a constraint minimization problem where we minimize the energy 
functional of the equation restricted to the subspace of functions that satisfy the conservation of 
mass restrictions. Then, in the Ritz sense, the obtained solution is the best among all functions that 
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satisfy the mass conservation restriction. 

As a main application of the techniques presented here, we consider the case where the coef¬ 
ficient k has high-variation and discontinuities (not necessarily aligned with the coarse grid). For 
this problem it is known that higher order approximation is needed. Indeed, in some cases ro¬ 
bust approximation properties, independent of the contrast, are required. See for instance lIsl-llOll 
where it is demonstrated that classical multiscale methods ( lIlSll l do not render robust approxi¬ 
mation properties in terms of the contrast. It is shown that one basis functions per coarse node 


9 , 


IJ. A similar issue 


1711 and related works. 


(with the usual support) is not enough to construct adequate coarse spaces 
can be expected for the multiscale finite volume method developed in [|15l - 
when applied to high-contras multiscale problems since the approximation spaces have similar 
approximation properties. In the case of Galerkin formulations, robust approximation properties 
are obtained by using the Generalized Multiscale Finite Element Method (GMsFEM) framework. 
The main goal of GMsFEMs is to construct coarse spaces for Multiscale Finite Element Methods 
(MsEEMs) that result in accurate coarse-scale solutions. This methodology was first developed in 
[3-3] based on some previous works [8-12]. A main ingredient in the construction is the use of an 
approximation of local eigenvectors (of carefully selected local eigenvalues problem) to construct 
the coarse spaces. Instead of using one coarse function per coarse node as in classical MsEEM, 
in the GMsFEM it was proposed to use several multiscale basis functions per coarse node. These 
basis functions represent important features of the solution within a coarse-grid block and they are 
computed using eigenvectors of an eigenvalue problem. Eor applications to high-contras problems, 
methodologies to keep small the dimension of the resulting coarse space where successfully pro¬ 
posed ([ 3 ])- That paper made use of coarse spaces that somehow incorporated important modes 
of a (local) energy related to the problem motivated the general version of the GMsEEM. The 


much more general GMsFEM was then developed in []5]] where several more practical options to 
compute important modes to be include in the coarse space were used; see also IllSn for an earlier 
construction. It is important to mention that the methodology in [3] was designed for parametric 
and nonlinear problems and can be applied in variety of applications, although a more extensive 
review of such developments is not contained herein. 

An important consideration of the proposed method is that the GMsFEM methodology does 
not guarantee conservation of mass properties such as those of Einite Volume formulations. In 
this paper we design a mass conservative GMsFEM method. In particular, we impose the mass 
conservation constraints over control volumes by using Eagrange multipliers. In doing so, we 
numerically validate that the convergence properties of the GMsFEM are maintained while the ap¬ 
proximate solution simultaneously satisfies the required conservation properties. We mention that 
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in [Il9n some successful applications of GMsFEM to two-phase problems were presented. In that 
work, the authors developed a technique in which GMsFEM solutions could be post-processed to 
yield mass conservative fluxes for use in two-phase modeling. While the motivation of the previous 
work in [Il9n is similar, the current technique yields a system that requires no post-processing, and 
a solution that automatically yields conservation. In particular, in the present work no additional 
equations must be solved (in notable constrast to lll9ll l in order to ensure coarse-grid conserva¬ 
tion, as a single global solve automatically yields the desired properties. As a result, no additional 
computational resources must be allocated, and issues such as ill-conditioning of the localized 
systems may be circumvented. Additionally, the proposed method ensures that the solution ob¬ 
tained through the method proposed in this work is the best among all functions (in the Ritz sense) 
satisfying the mass conservation restrictions. 

Finally we mention the works on MsFV. When dealing with multiscale problems, the lack of 
accuracy in the use of classical spaces had led researchers to introduce enriched spaces (as well 
as iterative procedures). In [|20n the authors construct an enrichment of the initial coarse space 
spanned by nodal basis functions. We also note that a hybrid finite-volume/Galerkin formulation 
for the coarse-scale problem is devised. The authors show that the resulting method has all fea¬ 
tures of the MsFV method, but is more robust and shows improved convergence properties if used 
in an iterative procedure. In comparison, in this work we use a GMsFEM framework were the ad¬ 
ditional degrees of freedom are related to local eigenvalue problems. This construction is related 
to the approximation properties the resulting space will have and is motivated by the numerical 
analysis of the interpolation operator into the coarse space. Furthermore, instead of a hybrid finite- 
volume/Galerkin formulation, we employ a Ritz formulation in the space of functions satisfying 
the conservation of mass. 

The rest of the paper is organized as follows. In Section [2l we recall how to impose restrictions 
using Lagrange multipliers in general and, in particular, for the mass conservation restrictions. In 
Section |3] we formulate a conservative coarse problem that follows the framework of Section [2l 
Section m is dedicated to present some relevant discussions for the case of smooth coefficients and 
solutions. In Section [5] we present the two-phase flow model problem and the overall solution 
algorithm. In Section we summarize some basic construction topics associated with the GMs¬ 
FEM methodology. In Section |7] we present some representative numerical results to illustrate 
the successful performance of our method for the case of single- and two-phase flow problems in 
high-contrast heterogeneous porous media. Some concluding remarks are finally offered in Section 
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2. Linear restrictions using Lagrange Multipliers 

Problem (O is equivalent to the minimization problem: Find p with (p — po) G Hjj and such 
that 


p = arg min J {v) 

V 


( 6 ) 


where the minimum is taken over v such that v — po E H\) and 


J{v) 


1 

2 


a{v,v) - F{v) + {gN,v)rj,. 


( 7 ) 


Let p be the solution of dH) and r*, i = 1,..., M, be M continuous linear functionals on H^{D). 
Define = Ti{p), i = 1, 2,..., M. The problem above is equivalent to: Find (p — po) G H]^ 
such that 

p = arg min J{v) (8) 


where 


W = {n : n — pd G and Ti{v) = rrii, i = 1,..., M}. 


Problem ([8]) above can be view as Lagrange multipliers min-max optimization problem. See 112 ill] 
and references therein. 

Then, in case an approximation of p, say p^, it is required to satisfy the constraints rj(p^) = rrii, 
i = 1,2,..., M, we can discretize directly the formulation ([8]). In particular, we can apply this 
approach to a set of mass conservation restrictions used in finite volume discretizations. 


2.1. Mass conservation in fine control volumes 

In order to ensure fine-scale conservation, we start by selecting control volumes We 

assume that each Vij is a subdomain of with polygonal boundary and i = l,...,M;.IfgGL2 
we have that dH) is equivalent to: Find p with {p — Pd) ^ Hjj and such that 


p = arg min J'[v) 
vew 

where the subset of functions that satisfy the mass conservation restrictions is defined by 


( 9 ) 


W = 


< n G H\) : f —AkVv ■ n= f q for all Vij > 

I -JaVij Jvij ’ J 


The Lagrange multiplier formulation of problem dH]) can be written as: Find p with {p — pn) G 
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Hjj and A G that solves, 


max min J7’(f ) — (a(p, /i) — F(/i)). (10) 

^eR^/ vGHj, 

Here, the average flux bilinear form a : H}j x —)■ M is defined by 

Mf 

= / -AkVv-n for all n G and/i G . (11) 

i=l JdVij 

The funetional F : —)■ M is defined by 


F(/x) = ^ /ij / q for all /i G 

The first order eonditions of the min-max problem above give the following saddle point problem: 
Find p with {p — pn) ^ Hjj and A G that solves, 

a{p, v) + a(r;, A) = F(v) - {qn, v)r^ for all v G 
a{p, p) = F{p) for all p G 


See for instanee 112 lh . 


3. Conservative discrete coarse problem 

Let Ff be a eoarse-mesh parameter and be a eoarse-seale triangulation. We assume that 
H does not neeessary resolve all the variation of the eoeffleient. In what follows we introduee 
the finite element spaee assoeiated with the eoarse resolution H. In applieations eoncerning 
heterogeneous multiseale media, the standard finite element spaees on do not offer good ap¬ 
proximation properties and therefore some speeial enriehment of this spaee is needed in order to 
aehieve some aeeeptable approximation properties. 

We denote by {yi}^Ai the vertiees of the eoarse mesh and define the neighborhood of the 
node Pi by 

G r^; y, G Kj}. (13) 

See Fig. [T]for an illustration of the eoarse-seale diseretization depleting eoarse neighborhoods and 
eontrol volumes. 

Using the eoarse mesh we introduee a set of eoarse basis funetions {$ 1 }. The basis fune- 
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Figure 1; Illustration of a coarse discretization depicting coarse neighborhoods and control volumes 

lions are supported in a;*; however, for one specific neighborhood uji, there may be multiple basis 
functions. We define the associated coarse space by 

Uo = span{$,}tT°^- (14) 

Let be a discrete interpolation of the boundary data pd- The classical multiscale solution of Q 
is given by pms with p^s — p% such that 

Pms = Siigmin J {v) ( 15 ) 

where the minimum is taken over the set of v with v — p^ G Vq. The equivalent matrix form of the 
Ritz-Galerkin formulation (fTSl) is given by 
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AoPms = b 

where the coarse matrix Aq and the coarse load vector b are such that 

u^Aqv = a{u,v) and v^b = F{v) — {gN, 


( 16 ) 

for all M, f e Uo we have 
v)t^. ( 17 ) 
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The multiscale solution p^s of problem (fT^ does not necessarily have any conservation of 
mass property. Therefore, the classical multiscale solution might not be appropriate form some ap¬ 
plications where it is imperative that approximations have some conservation properties in the form 
of average flux through some control volumes. To obtain multiscale solutions with the required 
conservation properties one can proceed as described in Section [2l To that end, we let Vq be de¬ 
fined as before and consider the set of all discrete functions that satisfy the required approximation 
properties. 

Wo = {v : V —p^ eVo and / —AkVv ■r]= q, i = 1,..., M^}, 

JdVi,, JVi,c 

where 1^ ^ denotes a coarse dual-grid volume. 

We then consider the following discrete formulation that takes into account the required re¬ 
strictions. The approximation of the solution of O is to find p/„ G Wo such that 

= arg min J"(n). (18) 

d€Wo 

This is a minimization problem with linear constraints similar to ([8]) - Indeed, it is the associ¬ 
ated Ritz-Galerkin formulation; see Section [2l This problem can be equivalently written in the 
form analogous to (fT^ or (fT^ . The matrix formulation stating the first order conditions of this 
minimization problem is written as 


Aq 

r' 


Pfv 


’ b ’ 

A 

0 


A 


b 


where the matrix Aq and the vector b are defined in (fTTI) and the finite volume matrix A is defined 
by 

fi^Av = a{v, /i) and fi^b = F{p) for all p G and v G Hl{D). 

4. The case of smooth solutions 

In this section, in order to motivate the use of the the formulation (flSI) for multiscale problems, 
we make some comments on the method described so far for the case of smooth coefficients and 
higher-order degree polynomials for finite volume methods. 

As mentioned in the introduction, when higher order degree polynomials are used within a 
finite volume framework one can proceed in different ways. Let us assume, only for this section 
and with motivation purposes, that Vq is the space of piecewise polynomials P''{T^) with r > 1. 
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Let us then introduee the dual grid T^*. The dimension of the approximation spaee Vq is higher 
that the number of eontrols volumes or dual elements. The number of eontrol volumes matehes 
the dimension of the subspace of piecewise linear functions P^(T^) contained in Vq. When r > 1 
and if we use only piecewise constant functions on the dual grid as a test function we do not obtain 
a square system. In order to obtain a square linear system one have to introduce additional test 
functions. One can proceed as follows. 


1. Construct additional control volumes and test the approximation spaces against piecewise 
constant functions over the total of control volumes (that include the dual grid element plus 
the additional control volumes). We mention that constructing additional control volumes is 
not an easy task and might be computationally expensive. We refer the interested reader to 
III 2. 22] for additional details. 


2. Use as additional the basis functions the basis functions that correspond to nodes other than 
vertices to obtain a FV/Galerkin formulation. This option has the advantage that no geo¬ 
metrical constructions have to be carried out. On the other hand, this formulation seems 
difficult to analyze. Also, some preliminary numerical tests suggest that the resulting linear 
system becomes unstable for higher order approximation spaces (especially for the case of 
high-contrast multiscale coefficients). 

3. Use the Ritz formulation with restrictions (flSl) . 


Note that if r = 1, in the linear system associated to (fT9l) . the restriction matrix corresponds to 
the usual finite volume matrix. This matrix is known to be invertible. In this case, the affine space 
Wo is a singleton. Moreover, the only function u satisfying the restriction is given by m = {A)~^b. 
The Ritz formulation (flSl) reduces to the classical finite volume method. 

If we switch to r > 1, the dimension of the space Wq can only increase. Therefore and due 
to ellipticity of the energy functional, the Ritz formulation with restrictions (ITSl) makes sense and 
it will give the best approximation of the solution in the space of functions that satisfy the restric¬ 
tions. Then, in the Ritz sense, the solution of (IT^ is not worse that any of the solutions obtained by 
the methods 1. or 2. mentioned above. Furthermore, the solution of the associated linear system 
(fT9l) . which is a saddle point linear system, can be readily implemented using efficient solvers for 
the matrix A (or efficient solvers for the classical finite volume matrix A); See for instance [|2in . 
Additionally, we mention that the analysis of the method can be carried out using usual tools for 
the analysis of restricted minimization of energy functionals and mixed finite element methods. 
The numerical analysis of our methodology is under current investigation and it will presented 


10 







elsewhere. 


The eomments and observations of this seetion are the main motivation for the methodology 
developed here to target fluid flow in porous media where, as it is well known, advaneed and 
sophistieated approximation spaees have to be used. Furthermore, in some applieations sueh as 
multiphase flow, eonservation of mass is a main requirement. 

5. Two-phase model problem 

We emphasize that solving the pressure equation in ([U) is done within the eontext of a two- 
phase flow model. In partieular, we are interested in treating a problem eonfined to a domain in 
whieh the subsurfaee is assumed to exhibit high-eontrast features. The heterogeneous reservoir is 
equipped with a well in whieh water is injeeted to displaee the trapped oil towards the produetion 
wells. The water and oil phases (whieh we denote by o and w, respeetively) are assumed to be 
immiseible, and we eonsider a gravity-free enviroment in whieh the pore spaee of the reservoir 
is fully saturated. Additionally, we assume that any eapillary effeets are negligible. Under sueh 
assumptions, eombining Darey’s law with a statement of eonservation of mass yields governing 
equations of the form 


V • V = g, 


( 20 ) 


where the total Darey veloeity is given by 


V = —A{S)k{x)'Vp. 

The pressure equations given by (l20l) is eoupled to a transport equation 

dS 


( 21 ) 


dt 


+ V-{vf{S)) = q^, 


( 22 ) 


where S denotes the water saturation, and q, q^ denote any external foreing. The total mobility 
ooeffleient A (S') and flux funetion /(S) are respeetively given by 

_ ^rwjS) ^ _ ^rw{S )//im 


Aw Ao 

where kra, a = te, o is the relative permeability of the a fluid phase. 
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Figure 2; An illustration of the operator splitting technique used to solve the two-phase model problem 


5.1. Solution algorithm 

In order to solve the coupled two-phase model given in Eqs. (I20l) . (|2TI) . and (l22l) we employ an 
operator splitting technique where the saturation from the previous time step is used to update the 
mobility coefficient required of the pressure solve and subsequent velocity calculation (see, e.g., 
11231]). Once this velocity is available, it is used in conjunction with an explicit saturation time¬ 
marching scheme for a specified number of time steps. The updated saturation is used again in 
order to update the pressure, and the process is continued until a final simulation time is reached. 
See Fig. [2]for a schematic of the the operator splitting. 

To discretize the saturation equation, we integrate Eq. (l22l) with respect to time, and then over 
some volume V) G O. Applying the left end-point quadrature rule to a second term in time, and 
integrating by parts yields the following expression 


meas(I/i)(S'i,„ - + At / v ■ n/(S'i,„_i) = / 

JdVi JVi 


where the errors terms have been neglected and 


5,- 


1 


S{x,tn). 


meas(I/i) Jy. 

In addition to the explicit time stepping in the above expression, we apply a non-oscillatory up- 
winding scheme on the flux term v ■ nf{Si^n-i) (see, e.g., ll24ll for a review of upwinding 
schemes on rectangular meshes). As mentioned in Section [2l it is crucial that that the numerical 
flux approximation satisfies a local conservation property. In particular, we wish to seek discrete 
solutions that satisfy 


/ Vft ■ n = 

JdVi 



for each Vi. 


(23) 


We emphasize that a main consideration of this paper is to ensure the conservation property in (|2^ 
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through the use of Lagrange multiplier restrietions (refer to Seetion ^ within a suitable eoarse- 
seale solution spaee (refer to Seetion [3]). In turn, in the next seetion we deseribe the systematie 
eonstruetion of the Generalized Multiseale Finite Element Method (GMsFEM) eoarse solution 
spaee. 


6. Generalized multiscale finite element method 


In this seetion we foeus on high-eontrast multiseale problems and summarize a GMsFEM 
eonstruetion of Vq used in Seetion |3l Eor a more detailed deseription of the development of the 
GMsFEM methodology, see 0 - 0 ] and referenees therein. 


6.1. Spectral enrichment 

We start by ehosing and initial set of basis funetions that form a partition of unity. The spaee 
generated by this basis funetions is enriehed using a loeal speetral problem. We use the multiseale 
basis funetions partition of unity with linear boundary eonditions (llSD). We have one funetion per 
eoarse-node and it is defined by 


— div(/cVxi) = 0 for iL e cuj (24) 

X* = Xi on OK, 

where is a standard linear partition of unity funetion. 

Eor eaeh eoarse node neighborhood ce*, eonsider the eigenvalue problem 

- div(A;) = afkijf , (25) 


with homogeneous Neumann boundary eondition on dui. Here and are eigenvalues and 
eigenveetors in oji and k is defined by 


k = kY,H^\^X3?- (26) 

i=i 

We use an aseending ordering on the eigenveetors, < aT < .... 

Using the partition of unity funetions from Eq. (l24l) and eigenfunetions from Eq. (1251) . we then 
eonstruet a set of enriehed multiseale basis funetions given by Xii^T seleeted eigenveetors 
Using Li to denote the number of basis funetions from the eoarse region Ui, we then define the 
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coarse GMsFEM space by 


Vo = span{<l>i,£ = i = i = 1,..., Li}. 


For more details, motivation of the construction, and approximation properties of the space Vq as 
well as the choice of the initial partition of unity basis functions we refer the interested reader to 

aa. 


Summarizing, in order to solve problem ([1]) for the pressure we use the GMsFEM coarse space 
Vo constructed in this section. Additionally, in order to obtain conservative solutions with respect 
to a dual coarse grid, we use the discretization presented in Section [3l In particular, we solve 
the saddle point problem (fT9l) using the appropriate approximation spaces and conservation con¬ 
straints. We also recall Section [5] for a description the overall solution algorithm used to solve the 
two-phase flow problem. 


6.2. A flux downscaling procedure 

At this point we carefully distinguish the scales at which we wish to calculate respective flux 
values. Recalling the discussions from Section 12.11 and Section |3l we note that sets of control 
volumes { or {Vi^c}ii\ are selected in order to incorporate the associated mass conservation 

restrictions of the fine and coarse problems, repsectively. In the case of the fine-scale problem, we 
simply assume that the set of control volumes {14,/} coincides with the nodal values of the mesh. 
In particular, the conservation constraints are taken directly on the fine scale. The case when we 
use GMsFEM is similar in the sense that we initially impose the constraints on the (larger) {I4,c} 
volumes associated with the coarse mesh discretization (see Fig. [U). As such, a global GMsFEM- 
FV solve directly yields flux parameters that satisfy a coarse analogue of conservation. While this 
level of conservation may hold for some target applications (see, e.g., [|25ll ). in this paper we wish 
to construct fluxes that are conservative directly on the fine scale for a direct means of comparison 
with the fine-sc ale FV technique. 

In constructing fine-scale mass conservative fluxes, we recall that the coarse solve from Section 
|3] yields the coarse conservation property 


Vff ■ n = 


q for all 14 ,c, 


JdVi,c 

which can analogously be thought of as the compatibility condition in 14,c- 
formulate a fully Neumann boundary value problem 


As such, we may 
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( 27 ) 


-V ■ [KkVpVi,,) =Q in Vi^c 
-KkVpVi,, ■ n = vh ■ n on dVi^c, 


where the known \^h = —AkVpfy is evaluated pointwise on the boundaries dVi^c- In particular, 
after obtaining the coarse solution and the coarse-scale conservative flux vh, we solve the set 
of localized problems in Eq. (l27l) for every 1/^ ^ in using any method that produces the desired 
fine-scale conservation. For consistency within this paper, we use the fine-scale FV technique to 
ensure fine conservation. A hallmark advantage of the downscaled post- processing procedure 
from (l27l) is that the localized problems are independent from one another. In particular, the prob¬ 
lems are naturally parallelizable (as should also be carefully noted for Eqs. (I24l) and (l25l)). and 
may be independently distributed to CPU and/or GPU multi-core clusters for a significant gain in 
computational efficiency [|26n . 


7. Numerical results 

In this section we offer a variety of numerical examples to test the performance of the method 
introduced in Section [2l In particular, we solve the model problem in Eq. ([T]) using the constrained 
discretization techniques from Sections [3] and [6l A main goal is to address the accuracy associ¬ 
ated with the reduced- order conservative GMsFEM approach as compared to the fine-scale FV 
approach. 

7.1. Single-phase pressure 

For the first set of examples, we employ both the fine-scale FV and GMsFEM-FV approaches 
to solve Eq. ([T]). Throughout the section we consider solutions that are obtained through solving 
the equation on the unit square domain D = [0,1] x [0,1]. We impose boundary conditions of 
Pl = I and = 0 on the left and right boundaries of the domain, along with no-flow (i.e., 
zero Neumann) conditions on the top and bottom. The fine coefficient (and reference solution) are 
posed on a 100 x 100 fine mesh that yields a global system of size Nf = 20402. See Fig. |3]for 
an illustration of high-constrast structure that is considered in this section. More specifically, the 
system from Eq. (fT9l) is of size Nf = dim(U^) -f Mf = 10201 -f 10201, where dim(U^) denotes 
the dimension of the space in which the fine-grid pressure solution is represented, and Mf denotes 
the number of dual-grid volumes where the finite volume constraints are imposed. In constructing 
the coarse-grid solutions in this section, we consider a 10 x 10 coarse mesh with varying levels of 
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Figure 3: A high-contrast permeability coefficient with inclusions and channels. 


constrained spectral enrichment. In particular, we obtain systems of size Nc = dim( Vq) +Mc where 
dim(Vb) denotes the number of degrees of freedom associated with the coarse enrichment, and Me 
denotes the number of coarse dual-grid volumes where the associated finite volume constraints 
are imposed. For this set of examples we emphasize that the number of coarse dual-grid volume 
constraints is fixed at Me = 121, which corresponds to the number of coarse nodal values (and 
surrounding volumes) of the domain (cf. Fig. [T]). 

As an initial motivation, we offer a comparison of single-phase pressure solutions in Fig. HI 


The coarse solutions in Figs. 4(b) and 4(c) were obtained through solving the global equation 


using two levels of coarse mesh enrichment. We can see from Fig. |4(b)| that a system of size 
Ne = 323 offers a somewhat crude approximation to the reference solution, yet that a system of 


size Ne = 647 (cf. Fig. |4(c)[ ) yields a solution that is nearly indistinguishable from the reference 
solution. In either case, we emphasize that the coarse systems are much smaller than the system of 
size Nf = 20402 that is used to obtain the fine-grid reference solution. 

For more rigorous comparisons we also consider the relative error quantities given by 


Ee = X 100%, (28) 

Ibll 

where p denotes the reference fine-grid solution, pe is a specified coarse-grid (GMsFEM-FV) so¬ 
lution. For the norm quantities in (|2^ . we consider both the energy and weighted norms 
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(a) Fine Nf = 20402 (b) Coarse Nc = 323 (c) Coarse Nc = 647 

Figure 4: A comparison between the fine pressure solution and increasingly accurate coarse-grid pressure solutions. 



(a) Fine Nf = 20402 (b) Coarse Nc = 323 (c) Coarse Nc = 647 

Figure 5: A comparison between the horizontal flux components of the fine system and the downscaled coarse system 


respectively given by 


\\p\\Hlin) = and \\p\\lUn) = . 

We note that the energy error is of particular importance for the analysis associated with GMsFEM 
(see, e.g., fl), and involves the gradient of the solution (which offers a flux- like comparison). 
However, using either norm we are primarily interested in illustrating the effects of larger coarse 
spaces (i.e., more basis functions) in the GMsFEM-FV construction. In Table[I]we offer a number 
of errors corresponding to a variety of coarse space dimensions. The left most column tabulates the 
total size of the full coarse system associated with a specified level of enrichment and conservation 
constraints. The next two columns itemize the dimension of the coarse space (and corresponding 
number of basis functions per coarse node), and the number of FV constraints. Most importantly, 
the results indicate that an increase in the coarse-space dimension yields a predictable error decline 
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Eull System 

Coarse Dimension 

Constraints 

Relative Errors (%) 

Nc 

dim(Eo) [# basis] 

Me 


mm 

242 

121 [1] 

121 

8.4 

>100 

323 

202 [2] 

121 

lA 

39.8 

485 

364 [4] 

121 

1.2 

15.0 

647 

526 [6] 

121 

0.9 

11.1 

809 

688 [8] 

121 

0.3 

9.0 

971 

850 [10] 

121 

0.3 

8.1 


Table 1: Relative pressure and flux errors for a variety of coarse space dimensions 


associated with the solution and its gradient. As sueh, these results strongly suggest that ineorpo- 
rating numerous pressure solves into the eontext of the operator splitting teehnique (reeall Seetion 
15.11) will yield inereasingly aeeurate two-phase saturation solutions. This is what we eonsider in 
the next subseetion. 

7.2. Two-phase saturation 

In this seetion we apply the eonstrained GMsFEM method to the full two-phase model as 
deseribed in Seetion [51 More speeifieally, we apply the eonservative GMsFEM diseretization of 
Eq. dl]) for eaeh pressure update. Then, a fine seale eonservative flux field is obtained via the 
downsealing proeedure in Seetion 16.21 in order to mareh the saturation solution in time using the 
explieit seheme from Seetion lSTl In doing so, we assess the effeetiveness of the proposed approaeh 
in a eontext where it is repeatedly used to aeeurately eapture a number of pressure equation updates. 

To solve the two-phase model given in Eqs. (l20l) and (l22l) we use quadratie relative permeability 
eurves given by krw = S'^ and kro = (1 — 5')^. We additionally use values of /i^ = 1 and /lo = 5 for 
the water and oil phase viseosities. The same pressure boundary eonditions from Seetion I tTI are 
used, and for the initial saturation eondition, we set S' = 1 at the left edge and assume S{x, 0) = 0 
elsewhere. Einally, we reeall that the high-eontrast permeability eoeffieient illustrated in Eig. [3] is 
used for k{x). Eor a motivating applieation of the method, we offer a representative set of two- 
phase flow solutions in Eig. [^ The plot shows saturation solutions advaneing in time for a variety 
of eoarse spaee dimensions, eompared with fine-seale referenee solutions. In partieular, the first 
row shows three saturation snapshots for the ease when Nf = 20402, and the seeond, third and 
fourth rows respeetively shows saturation snapshots for the eases when = 242, 323, 647. We 
note a signifieant improvement for the ease when Nc = 647 as eompared to the lower dimension 
Nc = 121. More speeifieally, the addition of more basis funetions to the eoarse pressure spaee 
yields flux and saturation values that aeeurately eapture the fine-seale behavior of the system. 
And, with respeet to the redueed dimension Nc = 647, we ean see that the solutions are nearly 
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indistinguishable from the referenee solutions. 

For a final set of eomparisons, we tabulate the standard relative error of the saturation 
profiles from Fig. as well as for a variety of other seenarios. The full set of results ean be seen in 
Fig. |7l For these examples, we use a variety of eoarse-spaee dimensions Nc = 242, 323,485, 647, 
and 809 and run the simulations to a final time of T = 0.9. A time stepping value of At = 10“"^ 
is used, and we update the saturation for 100 time steps in between eaeh pressure solve. As was 
evident from the previous illustration, we ean see from Fig. |7]that the error may be signifieantly 
deereased by adding more basis funetions in the GMsFEM eonstruetion. As a partieular example, 
a maximum error value of roughly 50% (Nc = 242) may be deereased to a maximum error value of 
roughly 8% in the ease when Nc = 809. These results serve to further illustrate the flexibility and 
aeeuraey assoeiated with the eonservative GMsFEM eonstruetion. More speeifieally, more basis 
funetions may be used to obtain a higher level of aeeuraey, whereas less basis funetions may be 
used when effieieney is a main eonsideration. 

8. Concluding remarks 

In this paper, we propose a method for the eonstruetion of loeally eonservative flux fields 
through a eonstrained variation of the Generalized Multiseale Finite Element Method (GMsFEM). 
The flux values are obtained through the use of a Ritz formulation in whieh we augment the result¬ 
ing linear system of the eontinuous Galerkin (CG) formulation in the higher-order GMsFEM ap¬ 
proximation spaee. We impose the finite volume-based restrietions through ineorporating a sealar 
Lagrange multiplier for eaeh mass eonservation eonstraint on a speeified seale, and as sueh, the 
proposed method may be viewed as a minimization problem in whieh an energy funetional of the 
governing equations is minimized within a subspaee of funetions that satisfy the desired eonser¬ 
vation properties. Due to the inherent eonstruetion of the eoarse-grid solution spaee, and the way 
in whieh the mass eonservation eonstraints are imposed, the eombined methodology is shown to 
offer a robust and flexible framework for obtaining eonservative flux fields to be used in two-phase 
flow modeling. To illustrate the performanee of the method we eonsider model flow equations 
with heterogeneous permeability eoeffleients that have high-variation and diseontinuities whieh 
signifieantly affeet the flow patterns of the two-phase model. The inerease in aeeuraey assoeiated 
with the eomputation of the GMsFEM pressure solutions is inherited by the flux fields and satura¬ 
tion solutions, and is elosely eorrelated to the size of the redueed-order systems. In partieular, the 
addition of more basis funetions to the enriehed multiseale spaee produees solutions that more ae- 
eurately eapture the behavior of the fine seale model. A variety of single- and two-phase numerieal 
examples are offered to validate the performanee of the method. 
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Figure 6; Water saturation profiles advancing in time for a variety of reduced-order dimensions 
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Figure 7; error results for the two-phase flow problem 
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